## Radon package
from numpy import *
from math import pi
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from FDMpkg import *
from scipy import interpolate
from optimalspec import *
from hilbert import *
def rickerwav(t,c,sig):
	""" A function generates a riker wavelet with a center at 'c' and
	sigma of 'sig'. The 't' should be a NumPy array containing the time grids. """
	y=2/(3*float(sig))**0.5/pi**0.25*(1.0-((t-c)/float(sig))**2)*exp(-0.5*((t-c)/float(sig))**2)
	return y
	
def gausswav(t,c,sig):
	""" A function generates a Gaussian wavelet with a center at 'c' and
	sigma of 'sig'. The 't' should be a NumPy array containing the time grids. """
	y=1.0/sig/(2.0*pi)**0.5*exp(-0.5*((t-c)/float(sig))**2.0)
	return y

def supergauss(t,c,sig,order):
	y=exp(-((t-c)/float(sig))**(2*order))
	return y	
	
def phaserot(x,phase):
	""" Rotate the phase by a constant phase of 'phase' of the inputed wavelet,
	'x' which is a NumPy array."""
	xf=fft.rfft(x)
	p=ones(xf.shape)*exp(1j*phase)
	p[0]=1
	l=p.size
	y=fft.irfft(xf*p)
	return y

def laplaceinterp(p,x,f):
	n=x.size
	re=0	
	for i in range(0,n):
		dn=1.0
		num=1.0
		for xx in x:
			if xx!=x[i]:
				dn=dn*(x[i]-xx)
				num=num*(p-xx)
		re+=num/dn*f[i]
	return re	

def shiftspec(x,shft):
	""" Time domain shift implemented by linear phase shift in spectrum domain."""
	xf=fft.rfft(x)
	l=xf.size
	ph=arange(0,l)*2*pi*float(-shft)/(2*l-1)*1j
	y=fft.irfft(xf*exp(ph))
	return y

def shift0pad(x,shft):
	y=shiftspec(x,shft)
	if shft>=0:
		y[0:int(shft)+1]=0
	else:
		y[int(shft)-1:]=0
	return y	
	
def profilechk(x,dot1,dot2):
	r=x.shape[0]
	c=x.shape[1]
	t=arange(0,r)
	if dot1[1]!=dot2[1]:
		a=(float(dot2[0])-dot1[0])/(dot2[1]-dot1[1])
		if a<=0:
			hmin=0
			hmax=min(c,ceil(dot1[1]-dot1[0]*a))
		else:
			hmax=c
			hmin=max(0,ceil(dot1[1]-dot1[0]*a))
		hrnge=arange(hmin,hmax)
		lc=a*(hrnge-dot1[1])+dot1[0]
		intpv=zeros(lc.shape)
		n=0
		for k in lc:
			s=max(0,floor(k)-3)
			e=min(r,floor(k)+4)
			intpv[n]=interp(k,t[s:e],x[s:e,hrnge[n]])
			n+=1
		plt.plot(hrnge,intpv)
		plt.show()
	else:
		print 'Two dots are at the same offset!'
	
def radonw(x,p,dh,dz,weight=None):
	""" Radon transform - 'x' is a 2D NumPy array containing the input 2D image, 
	'p' is a NumPy array containing p grids, 'dh' is the increment of rows of 
	'x', and 'dz' is the increment of columns of 'x'. """
	bd=5.0
	vsd=0.05
	derd=30
	y=zeros((x.shape[0],p.size))
	nh=x.shape[1]
	ch=(nh+1)/2
	tl=x.shape[0]
	t=arange(0,tl,1)
	if weight=='derivation':
		avamp=sum(abs(x))/x.size
		
	for tau in t:
		ip=0
		for pp in p:
			if pp!=0:
				if pp>0:
					hmn=max(ceil(-tau*dz/dh/pp)+ch,0)
					hmx=min(floor((tl-1-tau)*dz/dh/pp)+ch,nh)
				else:
					hmn=max(ceil((tl-1-tau)*dz/dh/pp)+ch,0)
					hmx=min(floor(-tau*dz/dh/pp+ch),nh)
			else:
				hmn=0
				hmx=nh
			count=arange(hmn,hmx,1)
			lc=tau+(count-ch)*pp*dh/dz
			intpv=zeros(lc.shape)
			k=0
			for c in lc:
				s=max(0,floor(c)-3)
				e=min(tl,floor(c)+4)
				intpv[k]=interp(c,t[s:e],x[s:e,count[k]])
				k+=1
			if weight=='boundary':
				if intpv[0]!=0:
					intpv[0]=intpv[0]*exp(-bd*abs(intpv[1]/intpv[0]-1))
				if intpv[-1]!=0:
					intpv[-1]=intpv[-1]*exp(-bd*abs(intpv[-2]/intpv[-1]-1))
			if weight=='derivation':
				deriv=diff(intpv)
				deriv=insert(deriv,[0],deriv[0])
				denom=abs(intpv)+avamp/1e8
				deriv=abs(deriv)/denom/(1+pp**2)**0.5
				intpv=intpv*exp(-derd*deriv)
			y[tau,ip]=sum(intpv)
			if weight=='variationsum':
				vsum=sum(abs(diff(intpv)))
				if tau==245 and pp==2:
					print abs(diff(intpv))
					print vsum
					print 2*abs(x[tau,ch])
				vsum=max(vsum-2*abs(x[tau,ch]),0)/(1+pp**2)**0.5
				if y[tau,ip]!=0:
					y[tau,ip]=y[tau,ip]*exp(-vsd*(hmx-hmn)*vsum/abs(y[tau,ip]))
			ip+=1
	return y

def radonm(x,p,hc,dh,dz):
    xf=fft.rfft(x,axis=0)
    N=xf.shape[0]
    L=xf.shape[1]
    J=p.size
    scl=zeros(xf.shape)
    R=ones((N,J))*1j
    for iw in range(0,L):
        scl[:,iw]=arange(0,N).T*2*pi/N
    ip=0
    for pp in p:
        ph=ones(xf.shape)*1j
        for iw in range(0,N):
            for ix in range(0,L):
                ph[iw,ix]=exp(1j*pi/(N-1)*(ix-hc)*iw*dh/dz*pp)
        R[:,ip]=sum(xf*ph*scl,axis=1)
        ip+=1
    y=fft.irfft(R,axis=0)
    return y

def iradonm(p,x,pc,dp,dz):
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=pf.shape[1]
    J=x.size
    R=ones((N,J))*1j
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*(pc-ip)*iw*dp/dz*xx)
        R[:,ix]=sum(pf*ph,axis=1)
        ix+=1
    y=fft.irfft(R,axis=0)
    return y

def radonslstk(x,p,hc,dh,dz):
	L=x.shape[1]
	N=x.shape[0]
	y=zeros((N,p.size))
	pmax=max(abs(p))
	rfm=ceil(dh*pmax/dz)
	dhm=dh/rfm
	ogrid=(arange(0,L)-hc)*dh
	fgrid=arange(-hc*dh,(L-1-hc)*dh+dhm,dhm)
	xc=zeros((N,fgrid.size))
	for i in range(0,N):
		intp=interpolate.interp1d(ogrid,x[i,:],kind='cubic',bounds_error=False,fill_value=0)
		xc[i,:]=intp(fgrid)
	ip=0
	xcc=zeros(xc.shape)
	for pp in p:
		for ix in range(0,fgrid.size):
			shft=-(ix-hc*rfm)*pp*dhm/dz
			xcc[:,ix]=shift0pad(xc[:,ix],shft)
		y[:,ip]=sum(xcc,axis=1)/rfm
		ip+=1
	return y

def iradonsprd(p,x,pc,dp,dz):
	L=p.shape[1]
	N=p.shape[0]
	y=zeros((p.shape[0],x.size))
	xmax=max(abs(x))
	rfm=ceil(dp*xmax/dz)
	dpm=dp/rfm
	ogrid=(arange(0,L)-pc)*dp
	fgrid=arange(-pc*dp,(L-1-pc)*dp+dpm,dpm)	
	pcp=zeros((N,fgrid.size))
	for i in range(0,N):
		intp=interpolate.interp1d(ogrid,p[i,:],kind='cubic',bounds_error=False,fill_value=0)
		pcp[i,:]=intp(fgrid)
	ix=0
	pcpc=zeros(pcp.shape)
	for xx in x:
		for ip in range(0,fgrid.size):
			shft=(ip-pc*rfm)*xx*dpm/dz
			pcpc[:,ip]=shift0pad(pcp[:,ip],shft)
		y[:,ix]=sum(pcpc,axis=1)/rfm
		ix+=1
	return y
	
def iradonhb(p,x,pc,dp,dz):
	y=iradonsprd(p,x,pc,dp,dz)
	L=y.shape[1]
	for i in range(0,L):
#		y[:,i]=-1/(2*pi)*dhilbertcheby(y[:,i],100,taper=0.1)
		y[:,i]=-1/(2*pi)*hilbertcheby(y[:,i],100,taper=0.1,coef=False)
	return y

def iradonmc(p,x,pc,dp,dz):
    dx=x[1]-x[0]		
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=pf.shape[1]
    J=x.size
    R=ones((N,J))*1j
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*(pc-ip)*iw*dp/dz*xx)
	    if iw!=0:
	    	prnge=floor((N-1)*dz/iw/dp/dx)		
	    	pmax=min(prnge+pc,L)
	    	pmin=max(pc-prnge,0)
	    else:
		pmax=L
		pmin=0
            temp=pf[iw,:]*ph[iw,:]
	    R[iw,ix]=sum(temp[pmin:pmax])
        ix+=1
    y=fft.irfft(R,axis=0)
    for ix in range(0,J):
	cnst=y[0,ix]	
	y[:,ix]-=cnst	
    return y


def radonm0(x,p,hc,dh,dz):
    xf=fft.rfft(x,axis=0)
    N=xf.shape[0]
    L=xf.shape[1]
    J=p.size
    R=ones((N,J))*1j
    ip=0
    for pp in p:
        ph=ones(xf.shape)*1j
        for iw in range(0,N):
            for ix in range(0,L):
                ph[iw,ix]=exp(1j*pi/(N-1)*(ix-hc)*iw*dh/dz*pp)
        R[:,ip]=sum(xf*ph,axis=1)
        ip+=1
    y=fft.irfft(R,axis=0)
    return y

def iradonm0(p,x,pc,dp,dz):
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=pf.shape[1]
    J=x.size
    scl=ones(pf.shape)*1j	
    R=ones((N,J))*1j
    for iw in range(0,L):
       scl[:,iw]=arange(0,N).T*2*pi/N
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*(pc-ip)*iw*dp/dz*xx)
        R[:,ix]=sum(pf*ph*scl,axis=1)
        ix+=1
    y=fft.irfft(R,axis=0)
    return y

def iradonm0c(p,x,pc,dp,dz):
    dx=x[1]-x[0]
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]	
    L=pf.shape[1]
    J=x.size
    pgrid=(arange(0,L)-pc+1)*dp
    R=zeros((N,J))*1j
    k0=zeros(J)*1j
    k0[J/2]=pf[0,pc-1]
    R[0,:]=roll(fft.ifft(roll(k0,-(J/2))),J/2)			
    for iw in range(1,N):
	prval=real(pf[iw,:])
	pival=imag(pf[iw,:])
	kgrid=-(arange(0,J)-J/2)*(2.0*(N-1)*dz/iw/J/dx)
	fr=interpolate.interp1d(pgrid,prval,kind='cubic',bounds_error=False,fill_value=0)
	fi=interpolate.interp1d(pgrid,pival,kind='cubic',bounds_error=False,fill_value=0)  
	kr=fr(kgrid)
	ki=fi(kgrid)
	R[iw,:]=roll(fft.ifft(roll(kr+1j*ki,-(J/2))),J/2)

    y=fft.irfft(R,axis=0)	
    return y

def iradonm0c1(p,x,pc,dp,dz):
    dx=x[1]-x[0]
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]	
    L=pf.shape[1]
    J=x.size
    pgrid=(arange(0,L)-pc+1)*dp
    R=zeros((N,J))*1j
#    k0=zeros(J)*1j
#    k0[J/2]=pf[0,pc-1]
#    R[0,:]=roll(fft.ifft(roll(k0,-(J/2))),J/2)
    wg=arange(0,J/2+1)				
    for iw in range(1,N):
	prval=real(pf[iw,:])
	pival=imag(pf[iw,:])
	kgrid=-(arange(0,J)-J/2)*(2.0*(N-1)*dz/iw/J/dx)
	dwg=wg*iw*dx*pi/(N-1)/dz
	if kgrid[0]<=pgrid[-1]:
		fr=fsinterp(pgrid,prval,dwg)
		fi=fsinterp(pgrid,pival,dwg)  
		kr=fr(kgrid)
		ki=fi(kgrid)
	else:
		fr=interpolate.interp1d(pgrid,prval,kind='quadratic',bounds_error=False,fill_value=0)
		fi=interpolate.interp1d(pgrid,pival,kind='quadratic',bounds_error=False,fill_value=0)  
		kr=fr(kgrid)
		ki=fi(kgrid)
	if isnan(max(kr)) or isnan(max(ki)):
		fr=interpolate.interp1d(pgrid,prval,kind='quadratic',bounds_error=False,fill_value=0)
		fi=interpolate.interp1d(pgrid,pival,kind='quadratic',bounds_error=False,fill_value=0)  
		kr=fr(kgrid)
		ki=fi(kgrid)
	R[iw,:]=roll(fft.ifft(roll(kr+1j*ki,-(J/2))),J/2)
#	print iw, max(kr)
    y=fft.irfft(R,axis=0)	
    return y
def iradonm0c2(p,x,pc,dp,dz):
    dx=x[1]-x[0]		
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=pf.shape[1]
    J=x.size
    R=ones((N,J))*1j
    scl=ones(pf.shape)*1j
    for iw in range(0,L):
       scl[:,iw]=arange(0,N).T*2*pi/N
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*(pc-ip)*iw*dp/dz*xx)
	    if iw!=0:
	    	prnge=floor((N-1)*dz/iw/dp/dx)		
	    	pmax=min(prnge+pc,L)
	    	pmin=max(pc-prnge,0)
	    else:
		pmax=L
		pmin=0
            temp=pf[iw,:]*ph[iw,:]*scl[iw,:]
	    R[iw,ix]=sum(temp[pmin:pmax])
        ix+=1
    y=fft.irfft(R,axis=0)
#    for ix in range(0,J):
#	cnst=y[0,ix]	
#	y[:,ix]-=cnst	
    return y
    	

def radonm1(x,p,hc,dh,dz):
    xf=fft.rfft(x,axis=0)
    N=xf.shape[0]
    L=xf.shape[1]
    J=p.size
    scl=ones(xf.shape)*1j
    R=ones((N,J))*1j
    for iw in range(0,L):
        scl[:,iw]=(arange(0,N).T*1j*2*pi/N)**0.5
    ip=0
    for pp in p:
        ph=ones(xf.shape)*1j
        for iw in range(0,N):
            for ix in range(0,L):
                ph[iw,ix]=exp(1j*pi/(N-1)*(ix-hc)*iw*dh/dz*pp)
        R[:,ip]=sum(xf*ph*scl,axis=1)
        ip+=1
    y=fft.irfft(R,axis=0)
    return y

def iradonm1(p,x,pc,dp,dz):
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=pf.shape[1]
    J=x.size
    scl=ones(pf.shape)*1j	
    R=ones((N,J))*1j
    for iw in range(0,L):
       scl[:,iw]=(-arange(0,N).T*1j*2*pi/N)**0.5
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*(pc-ip)*iw*dp/dz*xx)
        R[:,ix]=sum(pf*ph*scl,axis=1)
        ix+=1
    y=fft.irfft(R,axis=0)
    return y
	
def radonm2(x,p,hc,dh,dz):
    xf=fft.rfft(x,axis=0)
    N=xf.shape[0]
    L=xf.shape[1]
    J=p.size
    scl=ones(xf.shape)*1j
    R=ones((N,J))*1j
    for iw in range(0,L):
        scl[:,iw]=(-arange(0,N).T*1j*2*pi/N)**0.5
    ip=0
    for pp in p:
        ph=ones(xf.shape)*1j
        for iw in range(0,N):
            for ix in range(0,L):
                ph[iw,ix]=exp(1j*pi/(N-1)*(ix-hc)*iw*dh/dz*pp)
        R[:,ip]=sum(xf*ph*scl,axis=1)
        ip+=1
    y=fft.irfft(R,axis=0)
    return y

def iradonm2(p,x,pc,dp,dz):
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=pf.shape[1]
    J=x.size
    scl=ones(pf.shape)*1j	
    R=ones((N,J))*1j
    for iw in range(0,L):
       scl[:,iw]=(arange(0,N).T*1j*2*pi/N)**0.5
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*(pc-ip)*iw*dp/dz*xx)
        R[:,ix]=sum(pf*ph*scl,axis=1)
        ix+=1
    y=fft.irfft(R,axis=0)
    return y

def radonlsq(x,p,hc,dh,dz):
    N=x.shape[0]
    L=x.shape[1]
    J=p.size
    M=(N-1)/2
    R=zeros((L,J,2*M+1))
    for j in range(0,J):
        for l in range(0,L):
            t=(l-hc)*dh/dz*p[j]
            lc=M+t
            if ceil(lc)!=floor(lc):
                wf=ceil(lc)-lc
                wc=lc-floor(lc)
                flc=floor(lc)%(2*M+1)
                clc=ceil(lc)%(2*M+1)
                R[l,j,flc]=wf
                R[l,j,clc]=wc
            else:
                flc=floor(lc)%(2*M+1)
                R[l,j,flc]=1
    Rk=fft.rfft(R)
    xk=fft.rfft(x,axis=0)[0:M+1,:]
    qk=ones([J,M+1])*1j
    for k in range(0,2*M+1):
	yk=dot(xk[k,:],Rk[:,:,k])
        Hk=dot(Rk[:,:,k].conj().T,Rk[:,:,k])
        qk[:,k]=linalg.lstsq(Hk.T,yk.T)
    q=fft.irfft(qk)
    return q.T

def radonnw(x,xc=0,rs=8):
    xf=fft.rfft(x,axis=0)
#    savetxt("inputat1rowr.txt",real(xf[1,:]))
#    savetxt("inputat1rowi.txt",imag(xf[1,:]))	
    N=xf.shape[0]
    L=xf.shape[1]
    t=arange(0,L)-xc
    R=zeros(xf.shape)*1j
#    print R.shape
    p=arange(0,rs,2.0*rs/L)
#    print p.shape
#    print L
#    print p
    for k in range(0,N):
	if k==0:
		temp=sum(xf[0,:])/L
		R[0,:]=temp
	elif float(rs*k)/(N-1)>=1:
		R[k,:]=0	
	else:
		w=p*pi*k/(N-1)
		sp=optimspec(t,xf[k,:],w)
		Ls=sp.size
		if Ls==L:
			R[k,:]=sp
			print "The k is %f" % (k)
		else:
			R[k,:]=sp[1:Ls]
			print "The k is %f" % (k)	
#    plt.plot(abs(R[1,:]))
#    plt.show()
    y=fft.irfft(R,axis=0)
    return y

def radonnw2(x,xc=0,rs=8):
    xf=fft.rfft(x,axis=0)
#    savetxt("inputat1rowr.txt",real(xf[1,:]))
#    savetxt("inputat1rowi.txt",imag(xf[1,:]))	
    N=xf.shape[0]
    L=xf.shape[1]
    msp=abs(xf).max()
#   print msp
    t=arange(0,L)-xc
    R=zeros(xf.shape)*1j
#    print R.shape
    p=arange(0,rs,2.0*rs/L)
#    print p.shape
#    print L
#    print p
    for k in range(0,N):
	if k==0:
		temp=sum(xf[0,:])/L
		R[0,:]=temp
	elif abs(xf[k,:]).max()<1e-4*msp:
		R[k,:]=0	
	else:
		w=p*pi*k/(N-1)
		if float(rs)/L*k/(N-1)<1/float(L):
			sp=optimspec(t,xf[k,:],w)
		else:
			sp=optimspec4(t,xf[k,:],w)
		Ls=sp.size
		if Ls==L:
			R[k,:]=sp
			print "The k is %f" % (k)
		else:
			R[k,:]=sp[1:Ls]
			print "The k is %f" % (k)	
#    plt.plot(abs(R[1,:]))
#    plt.show()
    y=fft.irfft(R,axis=0)
    return y

def iradonnw(p,x,rs=8):
    pf=fft.rfft(p,axis=0)
    N=pf.shape[0]
    L=p.shape[1]
    J=x.size
    R=ones((N,J))*1j
    pc=L/2
    ix=0
    for xx in x:
        ph=ones(pf.shape)*1j
        for iw in range(0,N):
            for ip in range(0,L):
                ph[iw,ip]=exp(1j*pi/(N-1)*((ip-pc)*2.0*rs/L)*iw*xx)
        R[:,ix]=sum(pf*ph,axis=1)
	R[0,ix]=R[0,ix]/float(L)
        ix+=1
    y=fft.irfft(R,axis=0)
    return y
	
if __name__=="__main__":
	tau=250
	sig=8
	nh=40
	shft=2.5
	shft2=2
	shft3=2
	ratio=0.5
	hdrs=fdmhdr(0,159,25,0,1,25,0,512,4)
	h=arange(-nh+1,nh,1)
#	p=arange(-20,20.05,0.05)
	p=arange(-40,40.05,0.05)
	t=arange(0,512,1)
	tl=t.size
	h0=rickerwav(t,tau,sig)
#	h0s=shift0pad(h0, 300)
#	plt.plot(h0)
#	plt.plot(h0s,'r')
#	plt.show()
#	h0=gausswav(t,tau,sig)
	h00=rickerwav(t,tau-20,sig)
#	h00=gausswav(t,tau-40,sig)
#	h0=phaserot(h0,pi/2)
#	plt.plot(t,h0)
#	plt.show()
#	h00=phaserot(h00,pi/2)
	g=zeros((t.size,h.size))
	g[:,nh-1]=h0
	g2=zeros((t.size,h.size))
	g2[:,nh-1]=h0*ratio
	g3=zeros((t.size,h.size))
	g3[:,nh-1]=h00*ratio
	for k in range(1,nh,1):
		tmp=shiftspec(h0,shft*k)
		tmp2=shiftspec(h0,shft2*k)*ratio
		tmp3=roll(h00,shft3*k)*ratio
		if shft>=0:
			tmp[0:shft*k]=0
			tmp2[0:shft2*k]=0
			tmp3[0:shft3*k]=0
		else:
			tmp[tl-1-shft*k,tl]=0
			tmp2[tl-1-shft2*k,tl]=0
			tmp3[tl-1-shft3*k,tl]=0
		g[:,nh-1+k]=tmp
		g2[:,nh-1+k]=tmp2
		g3[:,nh-1+k]=tmp3
	for k in range(1-nh,0,1):
		tmp=shiftspec(h0,shft*k)
		tmp2=shiftspec(h0,shft2*k)*ratio
		tmp3=roll(h00,shft3*k)*ratio
		if shft>=0:
			tmp[tl-1+shft*k:tl]=0
			tmp2[tl-1+shft2*k:tl]=0
			tmp3[tl-1+shft3*k:tl]=0
		else:
			tmp[0:-shft*k]=0
			tmp2[0:-shft2*k]=0
			tmp3[0:-shft3*k]=0
		g[:,nh-1+k]=tmp
		g2[:,nh-1+k]=tmp2
		g3[:,nh-1+k]=tmp3
	gt=g
#	iq=gt
#	gth=genfromtxt("soff_gth1.txt")
#	print gth[:,0].size, h.size
#	gt=zeros((gth[:,0].size,h.size))
#	gt[:,nh-1:]=gth
#	gt[:,0:nh-1]=gth[:,nh-1:0:-1]
#	savetxt("inp_ricker_4hbt",gt)
#	fdmd=zeros((h.size,1,t.size))
#	fdmd[:,0,:]=gt.T
#	savefdm('radon_input_ricker.fdm',fdmd,hdrs)
#	plt.plot(abs(fft.rfft(gt[:,nh-1]))[0:200])
#	plt.show()
#	q=radonslstk(gt,p,nh-1,1,1)
#	x=arange(-80,81)
#	savetxt("ricker_40p_4hbt.txt",q)
        q=genfromtxt("ricker_40p_4hbt.txt")
#	q=q[:,7650:7851]
#	qc=q[:,39]
#	qfc=fft.rfft(q[:,29])
#	iq=iradonnw(gt,h,2)
#	savetxt("fake_gath.txt",iq)
#	iq=genfromtxt("inv_ricker_alias_gd.txt")
#	iqc=iq[:,39]
#	plt.plot(qc)
#	plt.plot(iqc,'--')
#	plt.xlabel('Time (a.u.)')
#	plt.ylabel('Amplitude (a.u.)')
#	plt.show()
#	iqfc=fft.rfft(iq[:,29])
#	iqf=fft.rfft(iq,axis=0)
#	savetxt("test_series.txt",iqf[50,:])
#	q=radonnw2(iq,79,rs=2)
	iq=iradonhb(q,h,800,0.05,1)
#	iq=genfromtxt("inv_ricker_40p_4hbt.txt")
#	savetxt("hinv_ricker_40p_4hbt.txt",iq)
#	iq=gt
#	J=qfc.size
#	dj=1.0/(J-1)
#	w=arange(0,1+dj,dj)
#	plt.subplot(211)
#	plt.plot(w,abs(qfc))
#	plt.plot(w,abs(iqfc),'--')
#	plt.ylabel('Amplitude (a.u.)')
#	plt.subplot(212)
#	plt.plot(w,angle(qfc))
#	plt.plot(w,angle(iqfc),'--')
#	plt.ylabel('Phase angle')
#	plt.xlabel('Frequency')
#	plt.show()
	vm=max(iq.max(),abs(iq.min()))
#	vm=max(q.max(),abs(gt.min()))
	print vm
#	plt.plot(abs(iqf[50,:]))
#	plt.show()
#	profilechk(q,[245,10],[230,30])
	plt.imshow(iq,cm.gray,aspect='auto',interpolation='bilinear',vmin=-7000,vmax=7000,extent=[-nh+1,nh-1,511,0])
#	plt.imshow(q,cm.gray,aspect='auto',interpolation='bilinear',vmin=-vm,vmax=vm,extent=[-8,8,511,0])
	plt.colorbar()
	plt.show()
